Computing Gibbs free energy differences by interface pinning 
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We propose an approach for computing the Gibbs free energy difference between phases of a 
material. The method is based on the determination of the average force acting on interfaces 
that separate the two phases of interest. This force, which depends on the Gibbs free energy 
difference between the phases, is computed by applying an external harmonic field that couples 
to a parameter which specifies the two phases. Validated first for the Lennard- Jones model, we 
demonstrate the fiexibility, efficiency and practical applicability of this approach by computing 
the melting temperatures of sodium, magnesium, aluminum and silicon at ambient pressure using 
density functional theory. Excellent agreement with experiment is found for all four elements, except 
for silicon, for which the melting temperature is, in agreement with previous simulations, seriously 
underestimated. 



An accurate location of first order transition lines at 
a reasonable computational cost is of paramount impor- 
tance for a wide spectrum of condensed matter systems, 
ranging from hard to soft materials and biological mat- 
ter. Basic principles of equilibrium thermodynamics im- 
ply that for a given temperature and pressure the system 
resides in the phase of lowest Gibbs free energy. Phase 
transitions occur where Gibbs free energy differences be- 
tween phases vanish, determining phase boundaries in 
the pressure-temperature plane. From the computational 
point of view, however, the task of evaluating a phase di- 
agram represents a significant challenge, as phase transi- 
tions occur on long time scales [T] such that they cannot 
be studied using straightforward molecular dynamics or 
Monte Carlo simulations. 

Several numerical approaches have been proposed to 
cope with this problem [H (i) in the indirect ap- 
proach, often based on thermodynamic integration, the 
Gibbs free energy is computed individually for each of 
the phases [IHS] and the coexistence line is then cal- 
culated by imposing the coexistence condition of equal 
Gibbs free energy, (ii) Alternatively, in the direct ap- 
proach, an explicit interface is introduced between the 
two phases which are then simulated simultaneously in 
the same simulation box. At fixed pressure and tem- 
perature, the system moves towards the phase with the 
lower Gibbs free energy. Exactly at coexistence the ther- 
modynamic driving force on the interface vanishes and 
the interface stops moving except for thermal fluctua- 
tions. Successful applications of this approach have been 
reported for a broad spectrum of materials [QtiST] . 

In this contribution, we present and validate a method 
to compute the Gibbs free energy difference, AG, be- 
tween two phases. The basic idea of this approach is 
to compute the average force required to pin the inter- 
face of a two-phase system via a harmonic bias potential. 
This external field couples to a suitably defined order 
parameter, Q, which distinguishes between the phases of 



interest. The application of the bias potential effectively 
transforms the out-of-equilibrium process of the conven- 
tional moving interface method into a well-defined equi- 
librium computation, in which the free energy difference 
AG is determined directly. We refer to this approach as 
the "interface pinning" method. Coexistence points may 
subsequently be determined using Newton's root finding 
method. 

To validate our new approach, we have first applied 
it to the Lennard- Jones (LJ) model[22]. Our calcula- 
tions reproduce with high accuracy the solid-liquid coex- 
istence line identified previously with other approaches 
[ini Hi] and provide Gibbs free energies that are in 
excellent agreement with those obtained by thermody- 
namic integration. We have then used interface pinning 
in combination with ah initio simulations to compute the 
melting temperatures of sodium (Na), magnesium (Mg), 
aluminum (Al) and silicon (Si), demonstrating that this 
method is efficient, flexible and widely applicable. 

Compared to the conventional direct and indirect 
methods used in the literature so far, interface pinning 
offers several advantages. In contrast to the direct ap- 
proaches, interface pinning operates at well-defined equi- 
librium conditions, thus permitting the explicit calcula- 
tion of free energy differences and interface properties. 
The selection of the order parameter Q does not need 
to be a reaction coordinate capturing the entire trans- 
formation mechanism. Finally, interface pinning inherits 
the general applicability and conceptual simplicity of the 
direct approaches. The latter makes it easy to implemen- 
tation into existing programs. 

To introduce the method, consider a two phase crystal- 
liquid system [25j in a periodic orthorhonibic box (see fig- 
ure [T]) at temperature T and pressure p. The box lengths 
X and Y are kept constant at values for which the crys- 
tal is unstrained, while the box length Z is allowed to 
change in order to maintain constant pressure. We refer 
to this ensemble as the TVp^T-ensemble. To lower the 
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interface Gibbs free energy Gi, the system will have two 
interfaces in the Xy-plane minimizing the interface sur- 
face area. We assume that the system is large enough 
to represent bulk phases at least at the center of the liq- 
uid and solid slabs. Particles may then either be labeled 
crystalline (subscript c), liquid (subscript I) or interfa- 
cial (subscript i), so that the total number of particles 
is N = Nc + Ni + Ni. The contributions to the total 
Gibbs free energy of particles in the bulk phases is deter- 
mined by the chemical potentials fic and /z/ of the solid 
and liquid, respectively, and the total Gibbs free energy 
is G = Nc^ic + Ni^ii + G,. 

When the relative distance between the interfaces 
changes, particles are transferred between the bulk 
phases. Assuming that the interface quantities Gi and Ni 
do not change when the interfaces shift due to the growth 
of one bulk phase at the cost of the other, the number of 
liquid particles may be written as Ni = —N^ + [constant] 
and the Gibbs free energy is given by 

G{Nc) = N^Afi + [constant] (1) 

where Afj, = Throughout the paper we will use 

the subscripts c and / for crystal and liquid properties, 
respectively, and let "A" denote "[crystal] — [liquid]". 

To sample configurations in the two-phase region and 
to prevent the system from complete transformation into 
one of the pure phases, we apply a harmonic bias poten- 
tial that pins the relative position of the interfaces. Let 
f7(R) be the energy of the unbiased system for configu- 
ration R = {ri, r2, . . . , rjv}, and 

C/'(R)=.;7(R) + |[Q(R)-a]2, (2) 

the energy of the system plus the bias potential. Here, 
(3(R) is a global order parameter with a linear depen- 
dence on the number of particles in the solid phase: 

Q = ^Qc + "j^Qi + 'kQ^ so that 

Nc = + [constant]. (3) 

In the biased system, the position of the interfaces rel- 
ative to each other will fluctuate around an average 
value and the order parameter Q will fluctuate accord- 
ingly. The probability distribution of Q is P'{Q) — 
exp[-G"(Q)/fcBr]/Z' where G'(Q) is the Gibbs free en- 
ergy along the Q coordinate of the biased system, and 
Z' is the corresponding partition function. The Gibbs 
free energy of the biased system may be written in terms 
of the unbiased free energy as G'{Q) = G{Q) + f (Q — 
af + kBT\n{Z'/Z). By insertion of Equs. ([l]) and ([s]), 
it follows that P'{Q) is Gaussian, 

where a = iVA/i/AQ is the slope of G{Q) in the two- 
phase region, displayed in the lower panel of Fig. [l] The 
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FIG. 1. (color online). Upper panel: crystal-liquid configura- 
tion from an ab initio simulation of 432 Si atoms in a periodic 
box. Atoms are colored according to the coordination num- 
ber (red= [fourfold coordinated] and blue otherwise). Lower 
panel: schematic sketch of the Gibbs free energy G{Q) (black 
solid line) as a function of the crystallinity order parameter 
Q at a state point where the liquid is thermodynamically sta- 
ble and the crystal is metastable. The double arrows indicate 
the interface contribution d (red) and the bulk contribution 
NAjj, (blue), respectively. The dashed green curve indicates 
the Gibbs free energy G'{Q) with bias potential applied. The 
inset shows that the computed G{Q) / {ksT) in the two-phase 
region is indeed linear; here, G{Q) was computed for the LJ 
model (A'^ = 5120) via umbrella sampling using Equ. 

S with K = 2 and a range of a's. 



distribution P'{Q) has variance ctq — kgT/K and average 
(Q)' = a — u/k, and the chemical potential difference 
between the two phases may be computed as 

Am = -kHQY - a)AQ/N. (5) 

As a guideline, we choose k, such that typical fluctuations 
in Q correspond to one or a fraction of a crystal plane, 
and a such that the system is approximately half liquid 
and half crystal. In practice, we find that a wide range 
of field parameters give the same precision of the A/i 
estimate [27] . 

Once A/i is known, coexistence points may be deter- 
mined using Newton's method for finding roots. The 
required derivatives of A/it along isobars and isotherms 
is given by the standard thermodynamic expressions, 
d{Aii)/dp\T = Av and d{A^i)/dT\p = -As = -{Au + 
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pAu — Aij,)/T. In these relations, Av, As, and Au are 
changes in specific volume, entropy, and energy, respec- 
tively. 

To apply the interface pinning method in practice, we 
must choose an order parameter Q that grows linearly 
with the number of crystalline particles Nc in the two- 
phase region. Moreover, Q should be computationally 
inexpensive. Unlike liquids, crystals have long-ranged 
translational order, allowing us to use the collective den- 
sity field as order parameter: Q = \pk\ where — 
X;^iexp(-ik-rj). Here, k = {2Tm^/X,27rny/Y,0) 
for some fixed integers (jixjUy) that should be chosen 
such that k correspond to a Bragg peak. This choice will 
maximize the contrast between the liquid and the crystal. 
The z-component of k is set to zero since Z fiuctuates in 
the A^pzT-ensemble. The constant N~2 makes Qi sys- 
tem size invariant while Qc oc . Derivatives of Q with 
respect to the particle coordinates, required to determine 
the forces resulting from the bias, can be computed with 
an algorithm scaling as 0{N). We note that this order 
parameter may be problematic in the supercooled regime, 
since a crystal can lower |pk| by introducing long wave 
length displacements of particles. The energy penalty of 
such displacements is low and decreases with increasing 
system size. We have chosen to use a-s order parame- 
ter for most computations, since it is generally applicable 
and simple. For some computations we have in addi- 
tion used the Steinhardt Q — Qq order parameter [5S], 
which has the advantages of being robust in the super- 
cooled regime. The two choices of order parameter give 
the same A/i's within statistical error. A more detailed 
description of the method will be given in a forthcoming 
publication j27| . 

To verify the method, we first used it to determine 
the solid-liquid coexistence line of the LJ model with 
truncated pair interactions: f/(R) = X]i>j- where 
u{r) = 4(r-i2 _ ^-6) _ 4(6-12 _ g-6) for r < 6 and 
zero otherwise (LJ units are used for this model through- 
out the paper). MD simulations with a time step of 
istep = 0.004 were performed using the LAMMPS soft- 
ware package [35] modified to include the bias potential. 
The Parrinello-Rahman barostat was used [33] with a 
time constant of tpr = 8 together with a Nose-Hoover 
[34), 135] thermostat with a time constant of tnh = 4. 

Results presented in Fig. [2] demonstrate that the solid- 
liquid coexistence line for the LJ model computed by in- 
terface pinning agrees within high precision with data 
obtained using other methods [MJ [30] . The coexistence 
points displayed in Fig. [2] were computed as follows. 
First, a crystal structure of 8x8x20 face centered cubic 
(fee) unit cells {N = 5120) was constructed and simu- 
lated at p = 1 and T = 0.8 for Un-a = 800. All box lengths 
were allowed to fluctuate in order to determine the ge- 
ometry of the unstrained crystal, giving X = Y — 12.85. 
The unstrained crystal was then simulated for tgim = 800 
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FIG. 2. (color online). Crystal-liquid coexistence of the LJ 
model (filled black dots) in the (p, T)-plane computed with 
interface pinning for the LJ model. The solid line is a cubic 
fit: -0.5223r^ + S.OITT^ + 5.502T - 5.989. The computed 
coexistence line agrees well with results of other methods |29) : 
-l-'s and x's are from Refs. [53] and |30| . respectively. The 
asterisk indicates the gas-liquid critical point (Tqp = 1.31; 
PQp = 0.15) of the full LJ model [31]. 



in the Np^T ensemble, and Qc = 56.31 {rix = 16, Uy — 0) 
and the average partial volume Vc = 1.036 was recorded. 
Next, a liquid was prepared by melting the crystal in a 
constant volume simulation at high temperature (T = 5). 
The iVp^r-ensemble (using X = Y ^ 12.85) of the hquid 
was simulated for tsim — 800, and Qi = 0.94 and the aver- 
age specific volume of the hquid vi = 1.163 was recorded. 
Then, a two phase configuration was constructed by per- 
forming a high temperature constant volume simulation 
where particles at z < Z/2 were kept at their crystal po- 
sitions using harmonic springs anchored at crystal sites, 
with the box volume (length Z) in between that of the 
crystal and the liquid. The iVp^T-ensemble with the 
bias-field of Equ. ([2| with parameters a = 26 and k = 4 
was simulated for igim — 4000 to compute (Q)' = 25.055. 
Application of Equ. ^ yielded a chemical potential dif- 
ference oi Aji — 0.040. The coexistence pressure was 
then determined iteratively using Newton's root finding 
method along the isotherm: — p(') — A/i*^*^/Ai;^'\ 

providing pressures of p^*^ = {1.0,1.320,1.337(1)}. In 
the last iteration the estimated chemical potential dif- 
ference is zero within the error bar, A/i = —0.0007(10) 
(throughout the paper numbers in parentheses indicate 
the statistical errors of the last digits) . In Fig. |3] we 
confirm that A/i(p, T) computed with interface pinning 
(symbols) is consistent with thermodynamic integration 
(lines). 

Due to its efficiency and fiexibility, the interface pin- 
ning approach can be combined with electronic struc- 
ture methods and ah initio molecular dynamics to com- 
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FIG. 3. (color online). Upper panels (a and c): Afi com- 
puted with interface pinning method along an isobar and an 
isotherm, respectively. Lower panels (b and d): specific en- 
tropy As vs. T and specific volume Av vs. p, respectively; 
the solid lines in the lower panels are quadratic polynomial 
fits to these data. The solid lines in the upper panels were 
computed by integration of these fits. The integration con- 
stants were chosen to provide the best overall agreement with 
the A/i-data. 



TABLE I. Ab mito and experimental (Refs. 1361 and [37 )) melt- 
ing temperatures Tm (in K) of period 3 elements using either 
\pk\ or Qe as order parameter. "Super cell" indicates the 
applied super cell built from the conventional cubic cell (in- 
cluding liquid and solid part). TV is the total particle number. 





unit cell 


super cell 


iV 


Q 


Tm [exp.] 


Na 


bcc 


5x5xfO 


500 


Pk & Qe 


354(21) [370] 


Mg 


hep 


4x6x8^ 


767 


Pk & Qe 


920(20) [923] 


Al 


fee 


4x4x8 


5f2 


Qe 


985(30) [933] 


Si 


cd 


3x4x7 


672 


IPkl 


1241(20) [1635] 



^ built from an orthorhombic 4 atom (a, \/3a, c) cell. 



puted free energy differences from first principles. For 
the present simulations, the method was implemented 
in the Vienna ab inito Simulation Package [38) . As an 
example, we used interface pinning to compute the melt- 
ing temperatures Tm of the period three elements Na, 
Mg, Al and Si at ambient pressure. Computed T^s are 
shown in Table |l] along with simulation details. Melt- 
ing temperatures were computed first for crystalline Si 
in the fourfold coordinated cubic diamond (cd) structure 
(see Fig. [T]). To be compatible to previous calculations 
im , density functional theory (DFT) in the local den- 
sity approximation (LDA) within the framework of the 
projector augmented wave method was used 00]. NpT 
and NpzT simulations were performed using a time step 
of tstep = 3 fs with a Parrinello-Rahman barostat [33] and 
a Langevin thermostat HTI. To compute the Si coexis- 



tence temperature at ambient pressure, we use a similar 
strategy as outlined for the LJ model: bulk properties 
of the crystal and the liquid {Qc, Qi, Vc, vi, X and Y) 
were evaluated in simulations for 216 Si atoms (3x3x3 
conventional cells; <sim = 60 ps) at T = 1200 K. Next, 
solid-liquid simulations with a bias field (tsim > 30 ps) 
were performed for four system sizes: {2x2x4, 2x2x7, 
3x3x6, 3x4x7} conventional cubic cells corresponding 
to iV = {128, 224, 432, 672} atoms. Coexistence tempera- 
tures were estimated to be {1189,1218,1225,1241} K us- 
ing Tm ~ T + Finally, finite size effects were ex- 
trapolated assuming a 1/N decay of the finite size er- 
ror yielding Tm = 1250(10) K. The finite size effects are 
particularly large for liquid silicon, since the metallic liq- 
uid is embedded in a semiconducting host, resulting in a 
discretization of the electronic states in the metal (elec- 
tron in a box) . The present value is fully consistent with 
previous LDA calculations [H |39], and the discrepancy 
to the experimental value of Tm = 1635 K originates 
from an underestimation of the energy difference between 
four fold coordinated semiconducting Si in the cd struc- 
ture and six fold coordinated metallic Si in liquid Si re- 
sembling the /3-tin structure |S]. The entropy of fusion 
As(Tm) = 3.5(1) /cs/atom and the slope of the melting 
curve dTra/dp — — 51(7)K/GPa (computed using the 
Clausius-Clapeyron relation) are also in agreement with 
previous theoretical results [HI [32] ■ For the other ele- 
ments, Na, Mg and Al, finite size effects are less critical, 
and we only considered system sizes comparable to the 
largest Si system. For these three elements, the calcu- 
lations were performed using PBEsol (Perdew, Burke, 
Ernzerhof functional for solids) ^42j . which yields more 
accurate lattice constants than the LDA. The computed 
Tm's of Na and Mg are in excellent agreement with ex- 
perimental values, while for Al the computed Tm is about 
6% too large (see Table |T]). 

In summary, we have introduced a computational 
method that allows a direct evaluation of the Gibbs free 
energy differences between two phases. In contrast to 
previous approaches, simulations are carried out at equi- 
librium conditions by pinning the interface between the 
phases of interest via a harmonic bias potential that cou- 
ples to a suitably defined order parameter. Application 
of interface pinning to the LJ model demonstrates the 
accuracy and efficiency of this new approach: the solid- 
liquid coexistence line agrees to high accuracy with data 
obtained by other methods and the computed Gibbs free 
energy is consistent with data obtained via thermody- 
namic integration. The practical applicability and flexi- 
bility of the method was demonstrated by computing the 
melting points of Na, Mg, Al and Si at ambient pressure 
using first principles simulations. The results demon- 
strate that present density functionals yield very accurate 
melting temperatures for crystalline metal to liquid metal 
transitions, but errors are sizable for semiconductor (cd- 
Si) to metal transitions (liquid Si). Furthermore, the 
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present approach allows to compute directly and straight- 
forwardly structural and thermodynamic properties of 
the interface, such as surface tension from the pressure 
tensor [?3] or the crystal growth rate from Q{t) fluctua- 
tions miis]. 
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